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Abstract 

We present a generalized discontinuous Galerkin method for a mul- 
ticomponent compressible barotropic Navier-Stokes system of equations. 
The system presented has a functional viscosity v which depends on the 
pressure p = p(p, in) of the flow, with the density p and the local con- 
centration pi. High order Runge-Kutta time discretization techniques 
are employed, and different methods of dealing with arbitrary coupled 
boundary conditions are discussed. Analysis of the energy consistency of 
the scheme is performed in addition to inspection of the relative error of 
the solution compared to exact analytic test cases. Finally several exam- 
ples, comparisons, generalizations and physical applications are presented. 

Keywords: Navier-Stokes; discontinuous Galerkin; Runge-Kutta; mix- 
ing; barotropic; compressible; viscous; miscible; multicomponent; multi- 
phase; multifluid; chemical; acoustic. 



Contents 

§1 Introduction 

§2 The generalized £-fluid 

§3 Towards a generalized boundary treatment 
54 Numerical test cases 



1 michoski@cm.ute.xas.edu, Department of Chemistry and Biochemistry 
* evans@ices.utexas.edu. Computational and Applied Mathematics 
1 pschmitz@math.utexas.edu, Department of Mathematics 
** vasseur@math.utexas.edu, Department of Mathematics 



1 



2 



Compressible Multifluid 



§5 Example: 2-fluid with chemical inlet 



§6 Example: fc-th order in time ^-fluid 



QH 



§7 Energy consistency of scheme 



§8 Fick's diffusion with acoustic BCs 



§9 Conclusion 



§ Acknowledgements 



§ Appendix 



1 Introduction 



Much work has been done in the study of the numerics of multicomponent 
flows. An example of an early yet comprehensive study of computational mul- 
tiphase mechanics was given by Harlow and Amsden in Ref. |25j . where they 
developed an implicit finite differencing technique for extremely generalized mul- 
ticomponent settings of both compressible and incompressible flows, including 
phenomena ranging from bubble formation and cavitation effects, to the for- 
mation of atmospheric precipitation and mixing jets. Subsequent and related 
work in multicomponent flows followed with, for example, the work of J. Dukow- 
icz in Ref. 18J for particle-fluid models of incompressible sprays, an approach 
extended by G. Faeth in Ref. [5U] to combustion flows and by D. Youngs in 
Ref. [55] to interfacial turbulent type flows. 

Owing to some of these pioneering works, recent work has demonstrated 
a resurgence of interest in multicomponent flows, approaches and numerical 
techniques. The importance of fluid-flows comprised of more than one phase, 
chemical constituent, species or component is represented by a vast array of 
applications that range across a number of fields. For example, multicompo- 
nent flows are essential for any flow demonstrating even rudimentary chemical 
kinetics; hence, for all (nontrivial) "chemical fluids" [55]. Likewise biological 
flows often require phase separations, in order to resolve membrane dynamics 
and interfacial behaviors in cells and cell organelles [15] and medical applica- 
tions desire estimates in local componentwise variations in blood serosity, which 
effect the viscosity and flow parameters involved with pulsatile hemodynamics 
[9 . Likewise we find numerous examples of multicomponent flow applications 
in the atmospheric [26] and geophysical (52] sciences; as well as in acoustics [39] 
and astrophysics [42] . just to mention a few. 

Here we present a new multicomponent numerical scheme based on a math- 
ematically well-posed [35] compressible barotropic system with functional vis- 
cosity depending on both the density p and the mass fraction \ii of each fluid 
component. It is well-known, both experimentally and theoretically, that viscos- 
ity has a functional relationship to the density and specie type (for examples see 
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the NIST thermophysical properties server). In addition, these types of mathe- 
matical models (with functional transport coefficients) are well understood from 
the point of view of continuum dynamics, having been extensively studied by 
Malek, Rajagopal et al. in Ref. [231 [32J [33J . It is further seen in Ref. [3H] that 
the analytic model used in this work a priori satisfies two essential entropy in- 
equalities, much like the shallow water equations [5], which serve as important 
tools for numerical analysis and implementation. 

In this paper we implement a discontinuous Galerkin (DG) finite clement 
method, employing piecewise polynomial approximations which do not enforce 
or require any type of continuity between the interfaces of "neighboring" el- 
ements. This particular implementation is primarily motivated by the works 
of Cockburn, Shu et al. (see Ref. [TTHTS] ) and Feistauer, Dolejsi et al. (see 
Ref. [16l[T7l[2lJ[22]). We implement a generalized formulation that is designed 
to accommodate an arbitrary choice of inviscid, viscous, and supplementary 
numerical fluxes. We use explicit time discretization methods as described in 
Ref. [12] . which necessitate a conditional stability requirement; namely the time 
discretization must satisfy the CFL condition. Up to the CFL stability con- 
dition we find our method to be very robust and to deal well with arbitrary 
numbers of fluid components of arbitrary type — up to the additional assump- 
tion that a barotropic pressure law is applicable. On the domain boundary data 
we again strive to generalize our setting. We show two different implementa- 
tions of boundary conditions, which demonstrate different solvency with respect 
to interior solutions, initial conditions and phenomenolgically relevant contexts. 
In both cases arbitrary Robin type BCs may be set. 

In §2 we give the general governing system of equations, the mathematical 
regularity, and the discrete formulation of the problem. In §3 we demonstrate 
a general way of dealing with boundary conditions by way of the method of 
characteristics, or alternatively, by way of setting arbitrary L°° data on the 
boundary. We provide an explicit formulation of the characteristic technique 
and show the generalized behavior of these types of "characteristic" boundary 
conditions, while subsequently discussing a number of alternative approaches. 
In §4 we implement two test cases with exact solutions, which are restrictions 
placed on the multifluid barotropic governing equations, and show that they 
are exact up to the possible exception of the boundary data. In §5 we show an 
example of a bifluid solution using the forward Eulcr method. We then show 
the difference between boundary conditions by way of weak entropy solutions 
versus that of characteristic boundary solutions. The next section, §6, is used 
to generalize the setting to n-fluid components and fc-th order Runge-Kutta 
schemes, where the example of an I — 5 fluid is shown explicitly. Then in 
§7 we analyse the energy consistency of the modelisation with respect to two 
entropy inequalities derived in Ref. |38j : one the classical entropy J? 9 and the 
second a closely related entropy 5^ discovered by Bresch and Desjardins (see 
Ref. [5, 6 ), where it turns out that the numerical scheme from §2 satisfies 
both energy relations provided the CFL condition is satisfied. Finally in §8 we 
extend the results to include Fick's diffusion law, where we inspect the exotic 
physical setting of a pressure wave traveling through a gas comprised partially 
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of polyynes, and discuss some applications. 



j2 The generalized ^-fluid 



We consider a one dimensional compressible barotropic ^-fluid system governed 
by the following system of equations: 

d t p + d x (pu) = 0, (2.1) 

d t (pu) + d x (pu 2 ) + d xP - d x {vd x u) = 0, (2.2) 

dtipfk) + d x {pupi) = 0, (2.3) 
with initial conditions, 

P\t=o = Pa > 0, pu\ t=0 = m , {ppi)\t=a = Pi,o- 

The multicomponent barotropic pressure p — p(pp*i, ■ . ■ ,pp n ) is chosen to sat- 
isfy 



P 



(ppiY 



(2.4) 



where J2i=i Pi — !• The mass conservation (2.11, momentum conservation 



(2.2 1, species conservation (2.3 1, and barotropic equation of state (2.4) describe 



the flow of a barotropic compressible viscous fluid defined for (t, x) € M + x M. 
Here the density is given as p, the velocity as u, the momentum by m, and the 
mass fraction of each component (chemical specie, phase clement, etc.) of the 
fluid is given by p,i, respectively, where > 1 corresponds to the emperically 
determined adiabatic exponent uniquely characterizing each of the £ species. 
Furthermore, adopting the notation throughout the paper that pi — ppi, the 
form of the viscosity functional v — v(pi, ■ ■ ■ , pi) is fixed to satisfy 



= ^'(Pi^PidpiP, 



(2.5) 



for ip'(p) = Cp a given a € (0, 1) and C > as emperically determined con- 
stants (see Ref. [2HE5]) and §7). 

The mathematical well-posedness of such a system (in the £ = 2 case) is 
given by the following theorem, which was proven by two of the authors in 
Ref. [38]: 



Theorem 2.1. Given (2.4) and (2.5) satisfying the conditions in Ref. I38f with 
initial data (pq,uq, (1q) satisfying 

< e(0) < Po < g(0) < ex), 
pa G H 1 (R), u GH 1 (R) , p,Q € i/ 1 (M), 

S'(po,po)dx < +oo, 
\d x p \ < Cp , 
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with ~g(0), g{0) positive constants and the i nter nal energy as given in Ref. 
there exists a global strong solution to (2.1)-(2.3) on K + x K such that for every 
T > we have 

perfO.T;^ 1 ^)), d tP e L 2 ((0,T) x R), 
u G L°°(0, T; H 1 (R)) n £ 2 (0, T; iJ 2 (M)), fyti G L 2 ((0,T) x E), 



Furthermore, there exist positive constants g(T) and g(T) depending only on T, 
such that 

< Q{T) < p(t, x) < q(T) < oo, V(t, x) G (0, T) x M. 

Additionally, when ip"(p), d pp p(p^), and d pfi p(p,[i) are each locally bounded 
then this solution is unique. 

Now notice that for an ^-fluid written with respect to conservation variables, 
the state vector U can be written as the transpose of the 1 x m row vector 

U = (p, pu, px, . . . p e ) T ', 

where m = I + 2 characterizes the degrees of freedom of our chosen system 
of equations. Note that we make this choice of a state vector for the sake of 
flexibility of representation and implementation (see for example §8), where the 
strict degrees of freedom of the system ( 2.1 1-( 2.3 ), due to the multiplicity of 
(2.1| in the conservation form of (2.3 1, is just I + 1. Nevertheless, consistent 
with our choice of an (£+2) x 1 state vector U, we obtain that the to x 1 inviscid 
flux vector / satisfies 



while the to x 1 viscous flux g is given by 

g(U,U x ) = {0,pu x ,0,. 



■,Peu) T , 

,of. 



In this notation (2.1 l-(2.3l can be expressed as 

U t + f x = g x , 



(2.6) 



where the notation (-) t corresponds to component- wise derivations with respect 
to the variable l. 

The Jacobian matrix of the inviscid flux Juf{U) = T(U) can be written as 
the to x to matrix: 



T(U) 






1 


••• 


-u 2 


2u 


d P iP ■ ■ ■ d Pt P 


-ti/Xl 


Pi 








ule 




Pn 





(2.7) 
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where li is the £ x £ identity matrix. An important feature of the barotropic 



pressure law ( 2.4 1 is that it is not a homogeneous function in pi, and thus the 
Jacobian T is not formulated to satisfy TU — f. This contrasts, for example, 
with the compressible Navier-Stokes equations when using the monofluid form 
of the ideal gas law p — Rpfl (see Ref. 22}). It should be noted that some 
numerical fluxes and schemes are designed or derived by specifically exploiting 
this homogeneity with respect to the Jacobian matrix of the flux function (for 
example, see the Vijayasundaram flux as used in Ref. [211 E2j). Nevertheless, 
our numerical fluxes will be defined independently of the homogeneity property 
of the corresponding map, where T simply satisfies f x = TU X . 



For the viscous flux g we define the m x m matrix 
Jf(U) = du x g(U,U x )=i 





_ u 1 
£ £_ 









1, (2.8) 







where here and below the O's are zero matrices of the appropriate sizes. Clearly 



then (2.6) satisfies 

U t +TU x -(J^U x ) x = 0. (2.9) 

Further let us introduce the auxilliary function X such that we are concerned 
with solving the coupled system: 

U t + TU X - (XY>) X = 

V ' 2.10 

s - u x = o. y ' 

The above equations comprise a first order system which can be effectively 
discretized using the DG method. 

Consider the following discretization scheme motivated by Ref. [22] (and 
illustrated in the one dimensional case in Figure [IJ . Take an open f2 C K with 
boundary Sf2 = T, given T > such that Qt = ((0,T) x fl) for f2 the closure 
of f2. Let denote the partition of the closure O, such that taking Ct = [a, b] 
provides the partition 

a — xq < x\ . . . < x ne — b 

comprised of elements Qi = (Xi—x,Xi) G S?h such that ^ = \Q\,Qi, ■ ■ ■ ,Gne\- 
The mesh diameter h is given by h = sxvpg e ar h (xi — Xi-i) such that a discrete 
approximation to Q is given by the set Qh = ^%Q% \ { a ,b}. Each element of 
the partition has a boundary set given by dQi = {ajj-i, x,-}, where elements 
sharing a boundary point dGi H dQj ^ are characterized as neighbors and 
generate the set JC-ij — dGi H dQj of interfaces between neighboring elements. 
The boundary d£l — {a, b} is characterized in the mesh as d£l = {xo,x ne } 
and indexed by elements Bj € d£l such that Cl = ,% L U K-ij U <9f2. Now for 
/ C Z + = {1, 2, . . .} define the indexing set r(i) — {j & I : Qj is a neighbor of 
Qi}, and for Is C Z~ = { — 1, —2, . . .} define s(i) = {j € Ib '■ Qi contains Bj}. 
Then for Si — r(i) U s(i), we have dQi — Uj e s(i)K>ij an d dQi H dfl = Uj Gs (j)/Cjj. 
We define the broken Sobolev space over the partition ^ as 

w k ' 2 (n h , sr h ) = {v : v lGz g w k ' 2 (Qi) vq,. e sr h }. 
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a b 
• • • • • • 



Xq X\ X2 X nG —2 X ne —\ X 



1 1 e 



o o o — • • ■ — o o o 

Ql Q2 Gne-l Que 



K12 /C 2 3 r fcne-l,ne 

^ne- 2,ne— 1 

Figure 1: The discretization of f2, distinguishing nodes, elements and neighbors, 
with boundary dfl — {a, b}. 



Further, approximate solutions to (2.1l-(2.3l will exist in the space of discon- 
tinuous piecewise polynomial functions over restricted to 3\, given as 

Sg(fi h> & h ) = {v : v lSi g <? d (G z ) VGi e ■%} 

for £? d {Qi) the space of degree < d polynomials on Qi. 

Choosing a degree d set of polynomial basis functions Ni G £? >d (Gi) for 
£ = Q, . . . , d we can denote an approximation to the state vector at the time t 
over Cl h , by 



U h (t,x)=J2ui(t)Ni(x), Vareft, 



where the Nfs are the finite element shape functions, and the U]'s correspond 
to the nodal unknowns. Likewise the test functions <p h: -&h G W 2,2 ^^, 37%) are 
characterized by 

d d 



1=0 



1=0 



for tp\ and i?; the nodal values of the test function in each Qi. 

Letting U be a classical solution to (2.101 and multiplying through by test 
functions ip h and and integrating clcmentwise by parts yields: 



d_ 
dt. 



\ U ■ <p h dx + (f ■ <p h ) x dx ~ f ■ (f^dx - g x - v h dx = 0, 

JQi JQi JQi JQi 

I T."d h dx- I (U-& h ) x dx+ I U-0%dx = O. 

JQi JQi JQi 



(2.11) 



Let <p\Kii an d PlJCji denote the values of ip on ICij considered from the interior 
and the exterior of Gi, respectivel. It should be noted that for ICij G V, the 
restricted functions <p h \ICji are determined up to a choice of boundary condition, 
which we will discuss in more detail in §3. Then we approximate the first term 
of pTT) by, 



d 

dt 



j g U h- <Ph dx ~ J t ' iPhdX ' 



(2.12) 
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the second term in (2.11 1 by the inviscid numerical flux 

*i(u h \K tj ,u h \ Kjt ,v> h )= 2^ / *(u h \ Ktj ,u h \ 



V / / • n i:j <p h \ Kij d1C, 



(2.13) 



for riij the unit outward pointing normal; and the third term on the left in (2.11 1 

by, 

&i(U h ,tp h ) = f f h ■ (<p h ) x dx « - / / • (<p h ) x dx. (2.14) 



The viscous term in (2.111 integrates by parts to give: 

/ 9x ■ fh dx = I (9- <Ph)zdx - / g ■ ip^dx 
JQi JQi J Si 



/ (JTS • (p h ) x dx - / JTS • ip^dx. 



(2.15) 



We approximate the first term on the right in (2.15 1 using a generalized viscous 
flux £f (see Ref. [T] for a review of choices in the DG framework) . We write here 
for the general viscous flux 

Kij^hWjiiU h\K i3 iU h\K oi ,nij) ■ ip h \/c i:j dlC 

• r nf -\ J Kin 



V / g ■ riijVhlKvdJC, 

^o/.s JKn 



(2.16) 



while the second term is approximated by: 

^ h ,U h ,ip h )= [ g h -(<p h ) x dxn [ g-ifldx. (2.17) 

JQi JQi 

Finally for the second equation in (2.10 1 we expand it such that the approx- 
imate solution satisfies: 



(2.18) 



J2 t (U,'E h ,U h ,# h ,#'>)= I V h -$ h dx + J U h -$ n x dx 



/ U{U h \ Kij ,U h \ Kji ,ti h \ Kij )dK, 



where, 



%{U h ,G h )= ]T I U{U h \ Kij ,U h \ Kji ^ h \ Kii )dK 
jes(i) K ^ 

f U- nij >& h \ Kij dK 



jes(i) JK >i 
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given that U is the generalized flux term associated with the discontinuous 
Galerkin method determined up to a congeries of options (please see Ref. [T] for 
a unified analysis), and using the approximate relations: 



X • -dhdx, and 



U h ■ ■d'ldx 



U ■ ■dZdx. 



We note that these choices of approximations and fluxes define the values of 
S/j on each clement in terms of the values of Uh on that element and adja- 
cent elements. As we shall see later, this indicates that with an explicit time- 
discretization scheme, computing is a completely local procedure. 
Combining \2.\ty , \2.\A\ , \2A&\ , ( |2.17| and ( |2.18| and setting, 



3C 



given the inner product 



i'lg 



E 



a,, 



bhdx, 



we define an approximate solution to (2.111 as functions Uh and for all 
t € (0, T) satisfying: 

l) E/fceC^M;^), X h eS(, 



2) j t (U h ,<P h )n g 



*{U h ,<p h )-G{U h ,<p h ) 

- 9(V h ,U h ,<p h ) + ^C£ h ,U h ,<p h ) = 0, ( 2 ' 19 ) 

3) £(tr,n h ,u h ,0 h ,0*) = a, 

4) U h (0) = U . 

We find below that up to a (possibly arbitrary) choice of boundary data, these 
solutions are quite well-behaved, extremely robust for arbitrary choice of n fluids 
(we show n = 1,2 and 5 here, and have tested up to n = 11 elsewhere) and 
readily extended to more complicated systems (e.g. §8). The results presented in 
this paper utilize piecewise linear basis functions, but we have tested quadratic 
basis functions in our code as well. 



§3 Towards a generalized boundary treatment 

Specifying arbitrary boundary data with respect to our approximate solution 



(2.19) is a delicate issue which requires a nuanced understanding of barotropic 
solutions and the mathematical techniques used to pose them. That is, we wish 
to determine the nature of an arbitrary boundary state U\qq in a way which is 



well-posed with respect to the system (2.1 )-(2.3); which is to say, in such a way 



that the uniqueness of the solution is maintained. 
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However, practically speaking, recovering boundary data of an arbitrary 
nature on d£l poses well-established difficulties with respect to the a priori es- 
timates established in Ref. [38], which serve as the cornerstone to the existence 
and uniqueness result stated in Theorem 2.1. That is, recovering the a priori 
estimates on the solution is reduced, in the first step, to recovering two entropy 
inequlities (see §7 for explicit forms) which serve as positive definite functionals 
over (0,T) x fl. However, when explicit boundary data is given, these inequal- 
ities acquire the addition of the following two unsigned boundary components, 
respectively (see Ref. [38] for the explicit calculation): 

(pu 3 ) x dx and / (pu(u + p" 1 ^) 2 )^, 
n Jn 

having the consequence of rendering the well-posedness of a formulation which 
spans any type of boundary data difficult to establish. Instead we offer a number 
of pragmatic approximate approaches that generalize the solution up to impor- 
tant restrictions, and then discuss some alternative approaches that are aimed 
at certain specialized types of settings. First we review some known results. 

It has been shown by Strikwerda in Ref. [50], and Gustaf 'sson and Sundstrom 
in Ref. [24 that incompletely parabolic systems, such as the shallow water equa- 
tions and the full Navier-Stokes equations, may be well-posed with respect to a 
broad set of initial-boundary data. These works additionally demonstrate the 
appropriate number of boundary conditions expected on incompletely parabolic 
systems, which differ from completely hyperbolic systems such as Euler's equa- 
tions. As the barotropic system ( |2.1| -(2.2 1 maintains a formal equivalence to the 



viscous shallow water equations (see for example Ref. [8, 37 ), we might expect 



(2.1 1- (2. 2 1 to behave as an incompletely parabolic system due to Ref. j23|. How- 
ever, the dependencies of the pressure p and viscosity v on the mass fractions 
make showing this nontrivial and require a careful analysis of either incompletely 
parabolic systems [50], or hyperbolic-parabolic systems [28] , 

The implementation of both incompletely parabolic and hyperbolic systems 
often rely upon the so-called "characteristic treatment." In these systems we 
use characteristic directions to extrapolate values of the system variables on 
the boundary, while the others become constrained by a set of characteristic 
relations (see Ref. [25] for the hyperbolic regime). These types of treatments 
have been extended to treat the full Navier-Stokes equations [13J 25J , the viscous 
shallow water equations [5U], and multifluid systems [ST] . 

We want to consider what we will refer to here and below as characteristic 
type boundary solutions, which we view as a reduced hyperbolic system (as 
presented in Ref. [22, 24.). We illustrate the situation for a simple one dimension 
case, but our analysis easily extends to the multidimensional case. To begin, 
suppose we have the domain (0, oo) in which to specify a characteristic boundary 
condition at the boundary point x = and time t = 0. Note that other one 
dimensional cases can easily be transformed to such a setting with a change 
of coordinates. We linearize our solution at x — with respect to a reference 
solution q, which for our purposes represents the numerical solution at x = 
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taken at a previous timestep. As an approximation, we neglect the viscous 
terms, resulting in: 

where q ~ U is a linearized approximation to the exact solution. Note that this 
arrives with a linear hyperbolic system. We consider the initial-boundary value 
problem in the set (0,oo) x (0,oo) equipped with the initial condition 

q(0,x) = q, for xe(0,oo), (3.2) 

and the boundary condition 

q(t,0) = q b (t). for te(0,oo), (3.3) 

Our goal is to choose the boundary condition q b (t) in such a way that the 
initial-boundary value problem is well-posed. To continue, we decompose into 
characteristic directions. That is, note that since T is diagonalizable we have 
that T/jCj = QCj, where the characteristic directions Cj are the eigenvectors of 
r associated to eigenvalues % (see §4 for an example). Then we can formulate 
the solution in the form 

m 

q{t,x) = J2 X j( t , x ) c j> ( 3 - 4 ) 

3=1 

where the initial and boundary data, respectively, satisfy 

m m 

Q = ^2 a i c 3> and q b = ^2^3 c r ( 3 - 5 ) 

3=1 3=1 



It then follows (from Ref. [52] chapter 3, for example) that (3.1) can be 
written as j initial-boundary value scalar problems: 

~dt + ^~d = {) in(0,oo)x(Q,oo), 

Xj(x, 0) = a 3 ;, for x £ (0, oo), (3-6) 

\ j {0,t)=(3 j fortG(0,oo). 



The scalar problems ( |3.6[ ) may be solved via the method of characteristics, from 
which we obtain the solution, 

^ ' for x - tq > 0, V ' 



which provides an explicit form to (3.4 1. From (3.7 1, we obtain the following 
conditions for the boundary data. 

• If > 0, then necessarily [3j — ctj. This is obtained by extrapolating the 
solution of Xj to the boundary x = 0. 
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Table 1: Choice of boundary conditions 



Boundary type 


Restrictions 


Free 


Fixed 


Subsonic 
u ■ n > — c 


inlet 
u ■ n < 


Mi H hMn = 1) P > 

and the appropriate 

suppiiiiiciiiai y 
conditions associated 
with a choice of 
boundary data, 
including: 
p,u,/j,i,p, u,m,pfj,i, etc. 




Pi 


Supersonic 
u ■ n < — c 


inlet 
u ■ n < 




none 


Subsonic 
u ■ n < c 


outlet 
u ■ n > 




01,03, 

■ ■ ■ j 0m 


Supersonic 
u ■ n > c 


outlet 
u ■ n > 


none 


01, 
■ ■ ■ , 0m 


Membrane 
u ■ n = 


wall 
u ■ n = 




01,03, 
, 0m 


Membrane 
u ■ n = 


osmotic 
u ■ n = 


02 , ■ ■ ■ , Pm 


01 



• If $j = 0, then Pj may be freely chosen. However, in some situations it 
may be useful to choose j3j = ctj for this case, such as an impermeable 
solid wall. 

• If $j < 0, then Pj may be freely chosen. 

Note that once we have selected well-posed characteristic boundary condi- 
tions, we utilize the transformation 

q b = V(0 1 ,P 2 ,0 3 ,--- ,p m ) T (3.8) 

to determine the consistent boundary conditions for the conservation variables. 



It turns out that for ( 2.1 1-( 2.5 1 we can reduce this method to that of the essential 
choices listed in Table [1] This corresponds with what we know of hyperbolic 
systems as shown in Ref. [55] and Ref. jlHji with respect to the number of free 
and fixed conditions on the boundaries. In Table [T] we also include a number 
of physically motivated restrictions which should be taken into account when 
selecting our boundary conditions. 

In addition to employing these "characteristic" solutions, we notice that the 



form of (3.1l satisfies the weak entropy solutions of Ref. [H|T5] and Ref. [33] for 
hyperbolic systems. However, as we show in §4, even though these two types of 
solutions are both consistent, they do not display equivalent numerical behavior. 
Nevertheless these two choices of boundary data, the characteristic and weak 



entropy solutions, are not ideal since (2.1 )-( 2.3 > is not a hyperbolic system 
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We may alternatively consider the route of positing boundary data by a simple 
extension of the results of Zlotnik (sec Ref. [58 ) to see that the barotropic system 
is parabolic in the sense of Petrovskii upon addition of the "quasihydrodynamic" 
or "quasigasdynamic" auxilliary function w (see Ref. [121157]) on <9D. 

More generally there exists a large back catalogue of results on compressible 
barotropic systems, many of which implement differing initial-boundary data, 
and some of which utilize fairly exotic conditions on the boundary. For example, 
for barotropic inflow problems we can refer to both Ref. [27] and Ref. (36] , where 
in both cases results from Ref. |38j are required and additional extensions are 
needed to move into the multiphase regime. Likewise solutions exist for free 
boundary barotropic problems jJH] , surface tension type boundaries [39] , Navier 
boundary type conditions |7J, and various Dirichlet type problems near vacuum 
states pHl EH 00] ; however, again, all of these results are only strictly satisfied 
for monofluidic systems, and thus require subtle analysis in order to extend 
them to the full multifluid regime. In many cases however, such as in Ref. |58j . 
the extension is relatively straightforward. 



§4 Numerical test cases 

We inspect two analytic test cases to verify the accuracy of the numerical method 
presented in §2 and §3. In both cases we solve a monofluid restriction of (2.1 1- 
(2.3 1 from the bifluid case (£ = 2), with p,\ = 1 and p,2 = in I = 1 spatial 
dimension. 

To begin, we specify the DG formulation in the bifluid case. First we define 
the three vectors U = (p, pu, p\, p>i) T , f(U ) = (pu,pu 2 + p, piu, p2u) T and 
g(U,U x ) = (0,vu x ,0) T such that (2.1 )-( 2.3 1 are expressed as 



u t + f x =g x 



whereby setting I = 2 in (2.7l and using ( 2 . 8 1 it then follows that 

u t + vu x = {,jeu x ) x . 



(4.1) 



(4.2) 



We can thus write a weak form of (2.1 1-( 2.3 1 in the same way as (2.11 1. 

To solve the system we must first specify the inviscid flux <I>. We test for 
two choices here. First we implement the local Lax-Friedrich's flux &ilf which 
satisfies 

/ $ilf ■ip h dK. = - / (f{U h )\ K +f(U h )\ K )-n ij <p h \ Kij dK 

-\f (Spec r (r))((t/ h )| K;i . - (U h ) Kji ) ■ ny^tadK, 

JKij 

for riij the outward unit normal and Spec r (r) the spectral radius of T. 
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As our second choice of inviscid flux we implement a standard approximate 
Ricmann solver, with flux & R satisfying: 



/ * R . Vh dK = \( (f(U h ) lKi] + (f(U h ) lfCji )-n %]l p h \ Kzj d)C 

m{Uh})\K{Uh})\V-\{U h })) ■ n^^dK, 



IK 

1 

~ 2 

where V and V~ l are found from the eigendecomposition given in the appendix, 
A is given by the diagonal matrix of eigenvalues diag(^) - as also enumerated 
in the appendix - and the average is given by 

{u h } = \{u h \^ + u h \ K]z ). 

Next we specify the viscous flux *S . Here we use a formulation similar to 
that presented in Ref. j4j , but we adapt it to include the functional dependencies 
present in the viscosity. We choose 

/ % ■ nijtp h dK = - I ((je"S h ) llCi) + {JfTYi h )\ Kji ) ■ u,,^ h K dJC. 

For the numerical flux U we use the Bassi-Rebay form, as shown in Ref. [TJ H] , 
which gives 

/ U BR (U h ,# h )dK = lf ((Uh)\K ti + (U h )\ Kji ) ■ ;<,.<>, , dk.. 

JKij 1 JKij 

Now we discretize in time, denoting a partition of [0,T] by 

= <° < t 1 . . . < t T = T, 

for a timestep given as At n = t n+1 — t n , and implement the forward Euler 
scheme: 

dU h _ u n h +1 - u n h 

dt ~ At" 

along with a slope limiting scheme in the conservation variables (p,pu), where 
the van Leer and Osher MUSCL schemes (as shown in Ref. |HJ [S3J |M]) have 
been adopted in this paper. 

Now we solve explicitly for (2.19). In particular, we show an explicit scheme 
using the Riemann flux, which is formulated to read: for every n > find U^ +1 
such that 

l)U%eS d h , K^si 



■jjn+l _ jjn 

2) ( h AM h ,<P h ) +* R (UZ,<p h )-®(U n h ,<p h , 



- u n h , <p h ) + u n h , Vh ) = o, (4 ' 3) 



3) 2(U BR ,2Z,Ut,0 h ,0Z)=O, 

h 




4) U h = U h (0) 



§4 Numerical test cases 
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Figure 2: The two graphs show the solution to (4.4 1 in terms of the linear Ric- 
mann flux and the Oshcr limiter, denoted u r and p r , versus the exact solution. 



The above formulation lends itself naturally to a staggered scheme. First, given 
Uh one solves step 3 for SJJ. This amounts to a simple, fast, and trivially 
parallelizable computation as the L 2 -projection matrix to be inverted is block- 
diagonal, with each block corresponding to an individual element. Second, given 
E^, one solves step 2 for U^ +1 . This similarly is a trivial computation as 
the mass matrix to be inverted is block-diagonal. In fact, with the choice of 
an orthogonal polynomial basis on each element, the L 2 -projection and mass 
matrices become diagonal. 

We inspect the first of two numerical test cases. Consider the monofluid 



steady state case of (|2.1|)-(|2.3|, by setting the initial data to po — fi x q = uq 



7, = 1 and p 2 ,o = 0- Clearly here the pressure reduces to unity p — 1 and the 
viscosity to a constant vq — Cq. Next we set the periodic boundary condition 

U n h {a+,t) = U n h (b+,t), U n h (a~,t) = Ut{b~,t). 

The exact solution shows constant solutions in the primitive variables. Our nu- 



merical simulations for 4.3 using both approximate Ricmann and Lax-Fricdrich's 
inviscid fluxes have shown that the L°° numerical error in the conservation vari- 
ables is of the order of machine precision, showing no fluctuation about the 
steady state in time. 

For the second of our test cases, we consider the monofluidic restriction of 



(2.1l-(2.3| given by taking /i^o = 1, p.2,0 = and = 1 with the additional 



relations: 

p = p = u^ 1 , and v = p. 
Solving this system immediately yields 

P^ 1 + P - pdxp^ 1 = -C, 
for Cg K, which leads to the ordinary differential equation 

u x = u 2 + 1 - Cu. (4.4) 
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Setting C = the noting that solution is independent of time, we solve the 
ODE yielding: u — tan 2. Setting the initial data to 

pa = (tanx) -1 , m = 1, and (p/ii)o = Pa, 

with the Dirichlet boundary data provided in the weak entropy sense of §3 via, 

p b = l/ub, m b = 1 and (ppi)b = Pb, 

we inspect the solution over the domain [a, b], with a = 0.5 and b = 0.7. Note 
that we enforce the weak entropy boundary conditions by setting U^k^ = 
{pb, m b , Pb, 0} in our discontinuous Galcrkin formulation. Here we compare the 
exact solution to the solution obtained using the Riemann flux with the Oshcr 
slope limiter (denoted p r in Figure [2J. 

In Figure [2] we plot the error of the numerical solution corresponding to a 
mesh size of h = 2 x 10~ 4 and a timestep size of h/30, where it is clear that the 
relative error over fifty timesteps is of the order of magnitude of the resolution 
of the mesh. The relative error is zero across the solution at the first timestep, 
as expected, and remains nearly constant in the interior of the domain in both 
cases, while the weak entropy implementation displays fluctuations in time of 
the order of h. These boundary fluctuations are neither monotonic nor generally 
increasing, but show complicated temporal perturbations at the weak entropy 
boundary points and are seen to weakly propagate into the interior as a function 
of the timestep. We have obtained similar behavior for the choices of a local 
Lax-Friedrich's inviscid flux and van Leer's slope limiter. Further, numerical 
experiments have revealed that the L 2 -error of the solution at a fixed time T 
scales like O(h) for the choice of a backward Euler scheme, a timestep size of 
At = h/30, and piecewise linear basis functions. For a general polynomial order 
d and an explicit time integration scheme of order k (see §6), we find the L 2 - 
error of the solution at a fixed time T scales like 0(h d+1 + At k ), as expected, 
provided the CFL condition is satisfied. 



§5 Example: 2-fluid with chemical inlet 

Let us show a simple application of the system outlined in §2 and §3 evaluated 
over two distinct constituents. Consider the bifluid system, 

d t p + d x (pu) =0, (5.1) 
d t (pu) + d x (pu 2 ) + d x p - d x (vd x u) = 0, (5.2) 
dt(pPi) + d x (pupi) = 0, (5.3) 

with initial conditions: 

p\t=o = Pa > 0, W|t=o = m and p, t=0 = fi . 

The pressure is given by p — pj 1 +pl 2 and the viscosity by v — ip' ("fipj 1 +72P2 2 ) 
for V' = Cp- a and a € (0, 1) with C > 0. 



§5 Example: 2-fluid with chemical inlet 



17 



Characteristic Inlet at (=12 




18 



27 



36 



• H,0 



45 



54 



Weak Entropy Inlet at / =12 
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Figure 3: The left plot shows miscible species at t = 12 given the characteristic 
chemical inlet conditions from (5.6 1 with = 0.9 and on the boundary a = 0, 
with the first order transmissive conditions on b = x ne (see Figure 0. The 
right plot shows the same solution using the weak entropy formulation. Here we 
have a miscible solution of methanol and water at $ = 500K and initial po = 5, 
u = 0, and /jig = fi2,o — 0-5. 



Now as in §4 we easily recover the form 

U t + TU x = (JtrU a ) s , (5.4) 



which integrates to (2.11 1. Again we solve for our system in a form equivalent to 
(2.19). We employ the local Lax-Fricdrich's inviscid flux $>ilf, the Bassi-Rebay 
numerical flux U br, the usual viscous flux and the van Leer slope limit er. 

All that remains is determining the boundary states U^\g u . We begin by 
considering the case of characteristic boundary conditions, and assume that at 
the boundary x = a we have a subsonic inlet u ■ n < 0. In our determination 
of characteristic boundary conditions, we linearize about the state Uh\Ku to 
arrive at an expression for the boundary state U^k^ at timestep t n . Then, 
from Table [TJ we see that /?" is fixed by 

f3[ l = p n (a+)/2, where p n (a+) = lim p n (a + x). (5.5) 

Now, suppose we want a chemical inlet such that the first chemical constituent 
fii is characterized by an influx condition /i™(a~) = c & where similarly, 

p n (a~)= lim p n (a — x). 
In order to maintain the consistency of our system, we additionally need that 



p^(a-) = \- c € and p n (a~) = e > 0. For I = 2, we solve (3.8 > with the 



constraint in (5.5) to obtain: 

ft = -p n (a + )/2-ft + e, % = ft? (o+ ) 



/3£ = eK(a+)-^)/4(a+), and (H"(0 = ft + ft + ft, '"^ 
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Difference in Characteristic and Weak Entropy 

0.08n 




-0.02 1 - 



Figure 4: Here we show the difference between the first chemical constituent 
of the weak entropy jjL\ tW and characteristic /ii jC solutions shown in Figure [3j 
where £ = fJ>i, w — fJ-i, c - I n this figure, for emphasis, we show only the reduced 
spatial interval (0,9). 

where £ n (a + ) = d Pl p n (a + ) — d P2 p n (a + ), c„(a + ) denotes the speed of sound at 
timestep t n on the boundary as defined in the appendix. Finally, at the other 
boundary point x = b we set a transmissive characteristic boundary condition. 

The behavior of such a "chemical inlet" is shown in Figure [3] where we have 
utilized a mesh size of h = 0.54 and a timestep of /i/30. Here we set e = po(a + ). 
By comparison the weak entropy solutions discussed in §3 to ( |5.1[ )-( [53| are also 
well-posed for an arbitrary collection of L°°((0, T) x dil) boundary data. So, in 
contrast to decomposing the solution into its characteristic directions, we may 
simply assign /z™(a~) = 'ta, ^(a^) = 1 — p n {a~) = e and the lag velocity 
condition u n (a~) — u n ^ 1 {a + ) for every timestep to obtain the weak entropy 
solution. 

Comparing the behavior of the weak entropy solution of the mass fraction of 
the first constituent in Figure [3] to the characteristic solution of the mass 
fraction of the first constituent \i\ tC yields Figure [4] Notice that the two bound- 
ary solutions do not demonstrate the same numerical behavior. In particular, 
the weak entropy y,\ grows more rapidly at the boundary; while the dynami- 
cally coupled characteristic solution adapts to the influx of specie/density by 
producing a velocity outflow, which effectively reduces the "chemical influx" as 
a function of time. 

In practice it is often physically meaningful to ascribe more boundary data 
than the free characteristic directions associated to the free /?'s can consistently 




§6 Example: fc-th order in time i'-fluid 
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Difference in Characterstic and Weak Entropy 




Figure 5: Here we have the complementary difference between species two of 
the weak entropy p,2,w an d characteristic pi,c solutions, where r\ = p,2,w ~ M2,c 

control. For example a closely related case to the chemical inlet example given 
above, is the subsonic outlet ■ n > 0. In cases such as these, where only one 
characteristic direction is free, weak entropy solutions are essential in order to 
even characterize such mixing at the boundary interface (such a case emerges of 
particular interest, for example, when interspecies diffusion occurs in the mass 
transport as shown in §8). Heuristically we may say that characteristic solutions 
demonstrate a relatively weaker forcing on <9f2 but are more restictive in terms of 
degrees of freedom, while the weak entropy boundary solutions display a greater 
flexibility of representation by way of establishing stronger forcing on dfl. 



§6 Example: k-th. order in time £-fluid 

We wish to generalize the example in §5 to £-fluid components and a k-th order 
in time Runge-Kutta time discretization. Let us start with an I = 5 system, 
which then can be easily generalized. Consider 

d t p + d x {pu) =0, (6.1) 
d t (pu) + d x (pu 2 ) + d x p - d x (vd x u) = 0, (6.2) 
dt(pni) + d x {pupi) = 0, (6.3) 

with initial conditions, 



P\t=o = Po > 0, pu\ t =o = m , and (pM»)|t=o = Pi,o, 
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Figure 6: Here we show the first and last timesteps of the mass fractions at 
•d = 293K using periodic boundary conditions with Runge-Kutta order k = 2. 
Initial conditions set p = 5+20e" (:l; " 10)5 V8+20e~( a!_30 ) 5 and u = sin(67ra;/a;„ e ), 
with fjti = 0.07 + 0.3e~ (K - 27 - 5)2/12 , fi 2 = 0.1 + O.Se^^- 10 5 ' 2 / 8 , p 3 = 0.06 + 
0.3e-( x ~ 22 - 5 '> 2 / 8 , fi 4 = 0.05 + o.Se-^- 30 - 5 ) 2 / 10 and solvent fi 5 = 1 - Yl=i Mi- 



given the pressure 

p = pT+pT + p? + pT + pT, (6-4) 

and viscosity 

v = ip'(pxd Pl p + p 2 d P2 p + p 3 d P3 p + P4d P4 p + P5d P5 p). (6.5) 

We take the three vectors U — (p, pu, pi, p 2 , p 3 , P4, p^) T , f(U) = (pu, pu 2 + 
p, piu, p 2 u, p 3 u, piU, p5u) T , and g(U, U x ) = (0, vu x , 0, 0, 0, 0, 0, 0) T , such that 
again we arrive with 

U t + TU X — (JffU x ) x = 0, (6.6) 



which is easily approximated by the numerical scheme given in (2.191. 

We generalize to higher order time discretization. That is, let us rewrite 

ferential equati 

U h = L h (U h ). 



(2.19 1 as a system of ordinary differential equations 

d 



dt 

We can solve this system using an explicit Runge-Kutta method. Specifi- 
cally, we use the strong-stability preserving Runge-Kutta methods presented in 
Ref. [H). This method follows for any €-fluid of the form ( |2.1[ )-(2.5 1 of Runge- 
Kutta order k. 

The behavior of this system is shown in Figures [6] and [7] where we have set 
the simple periodic boundary condition, 

U n h (a+,t) = U%(b + ,t), Ul(a-,t) = U n h {b~,t). 

The numerical solution shown was obtained using a mesh size of h = 0.36, a 
timestep of /i/30, the local Lax-Friedrich's flux, van Leer's slope limiter, and the 



§7 Energy consistency of scheme 



21 




Runge-Kutta method presented in Ref. [12] with k — 2. It is worth noting that 
the composition of this mixture does not tend towards homogeneous equilibrium, 
since there is both no interspecies diffusion (see §8) and the species are not 
"chemically miscible" (in that they do not mix in all proportions) . Nevertheless 
there is significant mixing from the state of the initial conditions, and it can 
be seen that the fluid is more homogenized, relatively speaking, at time t = 10 
that it was in the initial state. Most importantly, this scheme now immediately 
extends to an arbitrary ^-fluid. 



§7 Energy consistency of scheme 

In Ref. |38] it is shown that any solution for which the Theorem holds should 
satisfy two closely related entropy inequalities. The first, a classical integral 
inequality taking the form 

\t I {pu 2 + 2£}dx+ [ v\u x \ 2 dx < 0, (7.1) 

and the second owing to Bresch and Desjardins (see Ref. [5, 8J), as 

H I {p\u + p- 1 ^ x \ 2 + 2<?}dx+ I p-V^I^^O, (7.2) 

where the internal energy S = ^(pi, . . . , p n ) is specified as: 
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Entropy Comparison 
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Figure 8: Here we plot the integral forms S^t and S^t for C = 1 and a = 0.9, 
where J n M^dx and J n M^dx are represented by the first timestep. The spatial 
mesh is chosen with ne = 100 with At — 0.01. 



Entropy consistent numerical schemes are often formulated in the literature 



in order to explicitly enforce entropy inequalities such as (7.1 ) and (7.2 1 over all 
of Qt (viz. Ref. [3l|47]). For example enforcing (7.1l may be done by utilizing 



a change of variables of the conservation variable form of the state vector U, 
into the so-called entropy variable form W, which is achieved by writing the 
entropy functional Jff = pu 2 /2 + $ and then setting the state vector as the 
partial with respect to the conservation variables W = The difficulty of 

implementation of these energy schemes, which are inherently implicit methods, 
underscores the importance of conserving energy consistency of the solution, and 
further serves as motivation for testing how our explicit scheme behaves with 



respect to (7.1l and (7.2 1 



Here we inspect the entropy consistency of our scheme with respect to (7.1 1 



and (7.2 1 using the I — 5 fluid with periodic boundary data as shown in §6. 



From the numerical perspective, we expect our solution (2.191 to obey entropy 



consistency up to a restriction of the CFL stability condition, which for inviscid 
flows scale as Cih/Spec r T > At and for the complementary viscous flows like 
C2/1 2 / max(i/, 1) > At, where the CFL constants are characterized by C\, C2 € 
(0, 1). We note that we do not expect energy consistency for an arbitrary choice 
of boundary data. 



To examine whether the two inequalities (7.1l and (7.2| are satisfied, we 



first define ,%C = \p\u + p 1 d x ijj\ 2 + $ , and check that the spacetime integrated 



§7 Energy consistency of scheme 
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Figure 9: Here we compare the viscosity v, density p and pressure p of the 
periodic 5-fluid from §6 with a — 0.9, C = 0.5, 150 meshpoints and At = .006. 



functionals for our numerical solution satisfy: 



Mrdx 



\ 2 dxdt < I Jf Q dx, 
In 



(7.3) 



and 



y T = / Mrdx 



~ x i\}'\p x ^dxdt < I Jt dx. 



(7.4) 



We show the results of this calculation for an arbitrarily chosen set of param- 
eters in Figure [8] As is clear from the graph, both (7.3 1 and (7.4) are satisfied. 
In fact we have confirmed that (2.19) satisfies (7.3 1 and (7.4) whenever the CFL 



condition is satisfied, up to the choice of a constant. It is interesting to note that 
both of these inequalities are satisfied for an arbitrary choice of a and C in the 
numerical setting. This confirms that the mathematical result from Ref. [38J is 
substantially more restrictive than the numerical one. 

As a side remark, the functional behavior of the viscosity is a relatively 



unique property of our system (2.1 1-(2.3 1, which is to say that commonly com- 
pressible Navier-Stokes systems utilize constant viscosity coefficients (eg. see 
Ref. [22] chapter 4 and Figure [9| and thus the energy consistency and the CFL 
condition is not dynamically coupled to the solution components. However, for 
our system, since the viscosity is a function of time, the CFL condition must 
update to reflect the local viscosity magnitude at each timestep. 
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58 Fick's diffusion with acoustic BCs 



Although Theorem 2.1 only applies to systems of the form ( |2.1[ )-(2.3 1, the partic- 



ular numerical scheme outlined in §2 can be easily extended to more complicated 

systems; and indeed can be extended with similar numerical behaviors. As an 
example let us consider the 5-fluid, 

d t p + d x (pu) =0, (8.1) 

d t (pu) + d x {pu 2 ) + d xP - d x {vd x u) = 0, (8.2) 

dt{pni) + d x {pupi) - d x (p@id x p,i) = 0, (8.3) 

with initial conditions: 



P\t=o = Po > 0, pu\ t =Q = m , 



and 



(pMi)|t=o — Pi,o, 



given (6.4 1, (6.5 1 and ^ the diffusivity constants of each respective species. 
Here the system is equivalent to that in §6, except we have added the Fick's 
diffusion law term to the advection equation in p. Thus the state vector and 
inviscid flux remain unchanged, while the vector g becomes 

g(U, U x ) = (0, vu x , pS>ip x , p3> n pt) T , (8.4) 

such that the corresponding viscous flux matrix yields J(f — djj x g. 




Figure 10: A weak entropy solution to an oscillating pressure front propagating 
through a 5-component low density (~ 100 molecules per cm) gas at $ = 20K. 
The chemical constituents arc comprised of species found in dark interstel- 
lar molecular clouds, where representative fractional abundances are adopted 
and the solution space is appropriately scaled; with corresponding initial condi- 
tions: H 2 - 80%, He - 19.9%, and trace CO, H (atomic hydrogen), and HC 3 N 
(cyanoacetylene) . 



§8 Fick's diffusion with acoustic BCs 
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Figure 11: A characteristic solution to the same oscillating pressure front pre- 
sented inlTol 



We set an acoustic inlet condition, which is equivalent to identifying the 
sound pressure on <90. We suppose that the pressure on the boundary is a 
classical time-harmonic solution to the acoustic wave equation, namely, p b — 
Pa + ylosin(wi) for a driving amplitude Aq, an ambient reference pressure po = 
Si (/°oMi,o) 7i > an d an angular frequency U). 



Here we have solved (8.1 1- (8. 3 1 using a formulation which is meant to weakly 



mimic some of the conditions found in interstellar nurseries, or interstellar molec- 



ular clouds. The solution is shown in Figure 10 where it is notable that the 
traveling sound field p b dynamically responds to the changing speed of sound 
c throughout the medium - which scales like the root of the local change in 
pressure up to the local species concentration. The initial conditions and the 
diffusivities were estimated with the help of Ref. [35J HI] . 

For the case of weak entropy conditions, we set the lag velocity condition 
u b Wji ~ u b~ l \K-ij ano - determine the boundary values of the pi from their initial 
concentrations on cT2. Since in the barotropic case the total pressure satisfies 
Pb = pT j we then use the Newton-Raphson method to solve for roots in p^ 
of the following equation: 

I 

f(pb) = £(pi«,b) 7< ~ (po + A sinM)). (8.5) 

i 

This determines the values of p on the boundary, where Aq < po is the natu- 
ral positivity constraint on the pressure inlet. We allow antisymmetric inlets 
on dfl — {a, b} leading to the formation of supernodes within the fluid do- 
main. With our boundary data defined, we utilize the definition: U^Ik-i — 
{ Pb (t n ), m b (t n ), pi, b (t n ), p 2 ,b{t n ), P3, 6 (i n ), PA,b{t n ), P5,b{t n )}- The solution is plot- 



ted in Figure 10 for the domain (0,54), a mesh size of h = 0.135, a timestep 
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Figure 12: Here we plot the relative difference between the weak entropy pres- 
sure p w and the characteristic pressure p c . 

size of h/30, and the Runge-Kutta method of order k = 2. 

By comparison we solve the characteristic acoustic inlet boundary solution 
using the formalism presented in §3. As in §4 we linearize about the state 
Uh\K i: j to arrive at an expression for the boundary state U^Ku at timestep t n . 
Now, to determine well-posed characteristic boundary data we must dynamically 
switch between the five regimes (up to a choice of membrane condition for 
u ■ n — 0) listed in Table [l] since the pressure oscillation pulls the velocity 
between transonic inlet and outlet conditions. That is, we switch between the 
following cases: 

• Subsonic inlet: /3™ is fixed by V~ 1 q b , while (3%, ...,/?" are given by the 
equations fi^ b = /Lt i;0 , and X)i(P&Mi,fc) 7i - (Po + M sin(wt)) = 0, 

• Supersonic inlet: /3™, ...,/?" are given by the equations ^5 = u b = 



11 



b S and Ei(P&Mi,&) 7< ~ (Po + A)Sin(u;t)) = 0, 



Subsonic outlet: /3™, /JJ , . . . , /3™ are fixed by V 1 q , h , and we solve for 
by way of the pressure equation £\ {pb^i,b{02 )) 7i ~~ (Po + ^0 sin(wi)) = 0, 

• Supersonic outlet: /JJ , • ■ • , /3" are fixed by V _1 q b , 

• Wall: /3", . . . , /3™ are fixed by and we solve by way of the 
pressure equation Y^iiPbtM^iffii)) 7 ' - (Po + A sm(uit)) = 0, 

where we note that above we have set q b = U^\k ji . 

It can be confirmed by inspection of Figures [TO] and \TT\ that the characteristic 
solution demonstrates substantially sharper profiles than the analogous profiles 
in the weak entropy solution, and these peaks decay more rapidly in time. To 



show this more clearly, we display the difference graph in Figure 12 It is not 
clear a priori which solution is more phcnomcnologically predictive. 
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59 Conclusion 



We have shown an efficient and robust high-order numerical scheme for a mix- 
ing compressible barotropic viscous fluid comprised of up to I distinct chemical 
constituents. The DG solution was shown to be in very good agreement with 
two exact solutions derived by a choice of initial conditions, which demonstrate 
minimal numerical error at the weak entropy boundaries, as expected. The 
solution was then shown for two time-explicit schemes, the forward Euler and 
fc-th order explicit Rungc-Kutta schemes. Analysis of the method demonstrated 
the expected conditional stability up to a restriction by the CFL condition, and 
we further found that the numerical scheme up to this stability parameter is 
energy consistent, satisfying a novel entropy inequality; and that the energy 
consistency hold for a large family of physically relevant problems. We further 
provide a family of free boundary type solutions which are easily implemented, 
and which are numerically well-behaved, where either weak entropy or char- 
acteristic treatments are employed for comparative studies, and it is seen that 
indeed they demonstrate distinctly different behaviors even given (seemingly) 
equivalent initial data. 

A number of examples and potential physical applications were shown and 
cited in order to develop a sense of the large number of applications in chemistry, 
physics, engineering, and related fields. 

Future directions of the work include the expansion to higher spatial dimen- 
sions (2 and 3 dimensional meshes), the inclusion of Arrhenius type chemical 



equations to (2.3 1, the inclusion of temperature $ dependence into the model, 
the addition of fluid-structure interfaces, and the expansion of the modelisation 
to include ionic and polar species as well as dense plasmas (magnetohydrody- 
namic effects), surface tension and gravitational effects. 



Acknowledgements 



The first author would like to thank H. Gupta for providing references on ther- 
mophysical properties and additional insights into the interstellar media, and 
to further express sincere gratitude to Prof. J. F. Stanton for his continued 
support. The third author was partially supported by the Department of En- 
ergy Computational Science Graduate Fellowship, provided under grant number 
DE-FG02-97ER25308. The fourth author was partially supported by the NSF 
Grant DMS 0607953. 



28 



Compressible Multifluid 



Appendix 



We have that T is of the form 
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where we set for i = 1 , . . . , n the indeterminates Zj = d Pi p. Solving the char- 
acteristic equation det (T — Ik) = 0, the eigenvalues counted with multiplicity 
are, 

<Ti = u + c , <r 2 = u - c, £3 = u, q = u,... , <; n+2 = u 



where c = \J \i\Z\ + . . . + \i n Z n . While u has multiplicity n it is better to con- 
sider the eigenvalues in the three groups, u ± c, w, and the remaining (n — 1) 
copies of u as illustrated by the decomposition of the diagonalizing transforma- 
tion matrix 



V(U) = (c 1 ---c n ) 
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whose columns are the corresponding eigenvectors, which we abbreviate for 
convenience in the 3x3 block matrix form 
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where we have set X = (/X2, . . . , Hn) T and Y = (Z2, 
The inverse transformation matrix is given by 



■ j Z n ). 
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